Separation and fractionation of order and disorder in highly polydisperse systems 



O 

(N 



o 

C3 



I 

c 

O 

o 



(N 
> 

o> 
d 

On 
O 



X 



L. A. Fernandez, 1, 2 V. Martin-Mayor, 1,2 B. Seoane, 1,2 and P. Verrocchio 3, 2 

1 Departarnento de Fisica Teorica I, Universidad Complutense, Av. Complutense, 28040 Madrid, Spain. 
2 Institute de Biocomputacion y Fisica de Sistemas Complejos (BIFI), Spain. 
'" ' Dipartvmento di Fisica, Universitd di Trento, via Sommarive 14, 38050 Povo, Trento, Italy. 

We study a polydisperse soft-spheres model for colloids by means of microcanonical Monte Carlo 
simulations. We consider a polydispersity as high as 24%. Although solidification occurs, neither a 
crystal nor an amorphous state are thermodynamically stable. A finite size scaling analysis reveals 
that in the thermodynamic limit: a) the fluid-solid transition is rather a crystal-amorphous phase- 
separation, b) such phase-separation is preceded by the dynamic glass transition, and c) small and 
big particles arrange themselves in the two phases according to a complex pattern not predicted by 
any fractionation scenario. 

PACS numbers: 61.43.Fs, 62.10.+s,64.60.My 



I. INTRODUCTION 



Although in condensed matter physics spatial order is 
naturally linked to low temperatures, the presence of in- 
herently disordered interactions {quenched disorder) chal- 
lenges such scenario. The issue has been extensively ad- 
dressed in lattice systems (spin glasses, magnetic mate- 
rials in random field, etc..) where quenched disorder in 
fact inhibits spatially ordered structures (although not 
other types of order). Much less is known about off- 
lattice systems. The issue presents some practical con- 
sequences. For example crystallization of very viscous 
colloidal samples with size dispersion S, see Eq. (JTJ be- 
low, larger than 12% does not occur, even after several 
months spent from the sample preparation [lj . This leads 
to several basic questions about the equilibrium phase di- 
agram of polydisperse systems [3-fTol] . Does enough large 
polydispersity hinder crystallization? Is the glass phase 
stable rather than only metastable? Is there a dynamic 
interplay between crystallization and the glass transi- 
tion [3, |Tyj? And, probably at a more fundamental level, 
is thermodynamic equilibrium relevant at all to describe 
real polydisperse materials or these are instead inher- 
ently off-equilibrium over the experimental time scales? 
Answering such questions is crucial for condensed mat- 
ter physics, since polydispersity is found both in artificial 
(synthetic colloids, polymers) and natural systems, from 
supercooled liquids on the atomic scale up to biological 
fluids such as blood. 

An attempt to rationalize the experimental findings 
is the so-called terminal polydispersity scenario where 
a characteristic value St ~ 0.12 exists above which the 
homogeneous crystal becomes thermodynamically unsta- 
ble. There is not consensus however about what kind of 
structure should replace such single phase crystal. Den- 
sity functional analysis [7] predicts the instability of any 
crystal structure (even partial) above St, thus leaving 
the amorphous ones (either liquid or solid) as the only 
possibility. Yet, the moment free-energy approach [5[ 
predicts fractionation: phase separation between many 
crystal phases [though of the same ordering, face cen- 



tered cubic (FCC), for instance], each one with a much 
narrower size dispersion than S. Fractionation is sup- 
ported by a recent numerical simulation that found that 
a first-order fluid-solid transition actually occurs at any 
polydispersity [8j. However, the solid phase is quite com- 
plex, at least in the high polydispersity region. In fact, 
for S > 0.19 the transition regards only a fraction of the 
particles and the ordered state is inhomogeneous. Such 
state has been previously referred to as I-phaseQ. 

Here we study the high polydispersity region, in par- 
ticular the point S — 0.24. The corresponding S — /3 phase 
diagram (/3 is the inverse temperature, 1/T) is sketched 
in the inset in Fig. [3] This region is of great interest 
for various reasons. First, the amount of crystalline or- 
der for the coldest/densest configurations is unknown. It 
turns out to be phase-separated between a crystal and 
an amorphous state. The pattern of particle-size distri- 
bution among the two states does not follow any simple 
fractionation rule. Second, it has been suggested |8j that 
in this system the dynamic glass transition occurs in the 
stable rather than in the metastable fluid region. Our 
results support this claim in the large N limit. Besides, 
the detailed knowledge of the equilibrium structures is 
needed in order to design new experimental or numerical 
methods to drive the system towards such structures. 

The layout of the rest of this work is as follows. 
In Sect. |TT] we describe our model, the microcanonical 
ensemble (Sect. Ill Al) . and the considered observables 
(Sect. Ill Bp . Our simulation algorithm and our thermal- 
ization checks are described in Sect. IIIII Our main nu- 
merical results are described in Sect. IIVI We present our 
conclusions in Sect. |V] 



II. MODEL 

Take as a paradigm for polydisperse off-lattice sys- 
tems the polydisperse soft spheres (PSS) model. We 
consider particles of radius o~i , with i — 1,2,..., A. The 
particle size cr, is drawn from a probability distribution 
function (pdf) P(cr). Size polydispersity is in general 
characterized by a single parameter, S, defined as the ra- 
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tio among the standard deviation and the mean of Pier)'- 

(1) 



At least for small polydispersity, S seems to be the only 
feature of P(cr) that controls the physical results. 
Our particles interact via a pair potential: 



1 



V ij( r ) = e \ZT2+ X 
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X 12 -r 1212/13 



Vij (r) = if ajjj > x c , with a^i 



if Xij < x c , (2) 



, x c = 12T3. 



We take e as energy unit. The potential is basically the 
repulsive part of Lennard- Jones, 1/r 12 . The only role 
of the linear piece is to provide a smooth long distance 

cut-off [H G3. 

Our length unit, ctq, is fixed by 



al = J d(Xida j P(a i )P{a j ){(Xi + a j f 



(3) 



Although Q generalizes well known models for simple 
liquids 13] (one would then choose oo ~ 1 nm), the scale 
invariance of the 1/r 12 potential suggests that our model 
may describe as well colloids. For the colloidal case one 
would choose Co ~ 1 micrometer. In fact, the cutoff in 
the potential @ makes it short-ranged as it is appropri- 
ate for colloidal systems. 

We simulated N particles in a box with periodic 
boundary conditions at density p = cr^ 3 . Due to the 
scale invariance of the 1/r 12 potential, the thermody- 
namic parameter that controls the problem is the combi- 
nation r = pT- 1 / 4 (T is the temperature) j34j|. 

Here we study the case where the size distribution 
is flat (constant in the range [cr m i n , cr max ]). Sample-to- 
sample fluctuations are eliminated by picking the diam- 
eters in a deterministic way @, [l4| , 
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Observe that 

5 = 



1 (r-1) 
V3(r + lV 



with r = 



(4) 



(5) 



Hence, <7 max /er m i n — >• oo at Soo — l/\/3 « 0.57735. 

Since polydispersity hampers crystallization [l| , a glass 
transition is to be expected. Although most of this work 
has been performed in the microcanonical ensemble, let 
us mention that we have also estimated the glass tem- 
perature in the (N, V, T) ensemble by means of Monte 
Carlo (MC) simulations. We simulated the equilibrium 
fluid state using only standard Metropolis single-particle 
moves (different choices of microscopic dynamics lead to 
basically equivalent results, see 15]). To locate the ki- 
netic glass transition by computing the relaxation time r 



of the fluid for N = 500, 864 in the range T € [1.3, 1.46] 
(data not shown). Our definition of the kinetic glass tran- 
sition r g corresponds to the point when r surpasses the 
10 6 MC steps. Both for N = 500 and 864 particles, we 
find that T g = 1.455(5). 

The signification of r g is rather different, depending 
on whether one is studying liquids (i.e. <jq ~ 1 nm) or 
colloids (cto ~ 1 micrometer). In the colloidal 
standard MC step corresponds roughly to 0.01 seconds 
of experimental time [l6[, so that t ~ 10 6 MC steps ~ 3 
hours of physical time and r g corresponds to the experi- 
mental glass transition. On the other hand, for liquids 1 
MC step is roughly equivalent to one picosecond. Thus, 
t ~ 10 6 MC steps ~ 10~ 6 physical seconds, implying 
that . Tg ra ther corresponds to the Mode Coupling tran- 
sition [171 ]. Indeed, for most molecular and polymeric 
glass-forming liquids r at the Mode Coupling tempera- 
ture lies in the range 10~ 7 ' 5 and 10" 6 - 5 seconds fill- 



A. The constant energy ensemble 

We shall be working in the (N, V, E) ensemble. Specif- 
ically, we shall be using Lustig's microcanonical Monte 
Carlo [lj| in the formulation of [2p| . 

Let U be the total potential energy of our system, 



Thus, the total energy is 



E = U + K , (e = E/N) . 



(6) 



(7) 



where K = Yli^iPi ^ s the kinetic energy associated 
to the conjugated momenta {pj}. Here, we are consider- 
ing just one conjugated momentum per particle. As the 
kinetic energy is non-negative by definition, we should 
have E > U. The conjugated momenta are explictly 
integrated out (they are simply a conceptual device to 
introduce the ensemble (HI). 

A quantity of major importance in the microcanonical 
ensemble is the entropy density, sjv(e): 



exp[AsAr(e)] 



(27dV) 



N/2 



NT(N/2) 



(8) 



n 



i=i dn N_ 1 



9(e - u) 



The Heaviside step function, 9(e — u), enforces e > u. 
The microcanonical average of an arbitrary function of 
the particle positions {r}i and of the energy density e, 
0({r}i;e) is defined as 



(O)e 



/IL=i dr,0(M,;eK({r},;e) 
Xllili dr l w Ar ({r} 4 ;e) 



, (9) 



where. 



u> N ({rh;e) = (e - «)T-^(e - u) . (10) 
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B. Observables 



The inverse temperature 



The main observable in a microcanonical simulation is 
the inverse temperature, computed as a microcanonical 
expectation value at fixed energy e: 



(3(e) = 0)e, P = 



N -2 
2N{e - v 



Note that 



/3(e) 



dsjy(e) 
de 



(11) 



(12) 



The function /3(e) holds the key to connect the mi- 
crocanonical formalism with the canonical one. In- 
deed, the canonical probability density for e, Pp^ (e) oc 
exp[iV(sjv(e) — /3e)] can be recovered from /3(e): 



logP^ ) (e 2 )-logP^ , (e 1 ) 







N I de (/3(e) - /3) 

'ei 



(13) 

In the thermodynamically stable region (i.e. d/3(e)/de < 
0), there is a single root of /3(e) = /3, located at the value 
of e where '(e) is maximum. Instead, at phase co- 
existence there are several solutions for /3(e) = /3. Their 
interpretation is explained in Sect. IIV Al 



2. The particle- density field 



particles whose index i verifies \i — xN\ < 0.05N, hence 
only particles of similar size are considered): 



Qi(x) 



47T 

21 + 1 



i 

E 

m——l 



1/2 



\Qlm(x)\" 



(16) 



where (Yi m are the spherical harmonics): 



!lm{x) = 1CTF\> HmW 



E, <eJ(l )^(0 



iV 6 (i) 

^ ] ^im(^ij)' 



The index j in the latter sum runs over the Nf,(i) neigh- 
bors of the particle i and fij is the unit vector linking the 
position of particles i and j. Particles i and j are said 
to be neighbors if \\n — Tj\\ < A. In order to meaning- 
fully fix the scale A, we considered the average number 
of neighbors as a function of A, finding a plateau. The 
height of the plateau is remarkably /V-independent, but 
its width increases with N (so, the particular choice of 
A becomes less critical as N grows). We fixed the value 
A = 0.35 (in units of the maximum cut-off for the poten- 
tial 2er max x C ut), that lies in the plateau for all our values 
of N and for all our energies in the solid phase. 

Since we let the fraction of particles be a finite fraction 
x of N the Qi's are intensive quantities. In amorphous 
phases Qi(x) decrease as 1/y/N while in crystalline ones 
they remain of order 1. In particular, we consider the 
case I = 6. 



As we mentioned in the Introduction, we expect large 
particle-density fluctuations. In order to detect them, we 
study the Fourier-transformed density field at the small- 
est, non-vanishing wavenumber allowed by the periodic 
boundary conditions: 

T= \ (|p(2vr/L,0,0)| 2 + permutations) , (14) 

where L is the linear dimension of our cubic simulation 
box and the Fourier field is 



(15) 



Note that p(q), a function of the particles configura- 
tion, yields the static structure factor through S(q) — 
N(\p(q)\ 2 ) . In particular, p(0) is our non-fluctuating 
particle density p. 



3. Crystalline order parameters 

In order to study simultaneously crystallization and 
fractionation, we generalize the (rotationally-invariant) 
crystal order parameters (2li ,22] by measuring the crys- 
tal order only of a given set of particles X(x) (namely, 



III. NUMERICAL ALGORITHMS AND 
THERM ALIZATION TESTS 

In order to study the fluid-solid phase transition we 
implement a microcanonical MC strategy(l9l.[20|. Fixing 
the total energy density e, while the temperature and 
the potential energy fluctuate (see Eq . (fTTj) ) . we follow 
the evolution from one phase to the other by studying e 
in the energy gap between the two phases. This strategy 
turned out to be essential to assess the first-order nature 
of the phase transition in disordered Potts models [23]. 

The peculiarity of the polydisperse models addressed 
here, as compared with Potts and similar models, is 
in that the phase transition actually corresponds to 
a phase separation. In fact, our low energy state is 
inhomogeneousQ. Thus moving e from large values 
(fluid) to small ones (partly solid) we gently accompany 
the system during the growth of the spatially segregated 
regions. Although internal energy will not be the only 
reaction coordinate (see below), we have found useful 
to combine microcanonical MC with a modified Paral- 
lel Tempering (PT) algorithm [13, HI]. 

For the sake of clarity, we divide the remaining part 
of this Section in three paragraphs: particle move- 
ments at fixed energy (Sect. IIII A[) . Parallel Tempering 
(Sect. HiTB|) . and thermalization checks (Sect. ITlTC|) . 
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A. Particle movements at fixed energy 

The particle moves at fixed energy were, with 50% 
probability, either standard Metropolis single-particle 
moves, or global swap attempts (modified for a poly- 
disperse system). Let us recall that in a swap move, 
one attempts to exchange the position of two particles 
of different sizes (26[. Both for single-particle and for 
swap moves we compute the ratio of the microcanonical 
weights, defined in Eq . (flUl) . for the new and the old con- 
figuration u;^ d / ■ The new configuration is accepted 
with Metropolis probability min{l, w^ d /w^ ow }. 

To fully describe the swap algorithm, we need to dis- 
cuss how we choose the pair of particles, A and B, whose 
position we are trying to interchange. Note that one 
needs to balance two effects in polydisperse systems. The 
acceptance is larger the closer the two particle sizes are. 
However, exchanging very different particles produces a 
more significant effect when trying to equilibrate the sys- 
tem. Our compromise has been the following. We pick 
particle A with uniform probability over the N possibili- 
ties. We pick B with uniform probability among particles 
such that \<jb — &a\ < 0.2(<T max — <7 m j n ) . Particle B is ac- 
cepted with probability 1 if \gb — <Ja\ > 0.1(a max — cr m ; n ) 
or with probability 0.2 in the opposite case. In case of 
rejection, a new particle B is selected until a suitable 
candidate is picked. 

On the coexistence-line, swap moves reduce by three 
orders of magnitude the tunneling time between the fluid 
and the solid phase. 



B. The microcanonical parallel tempering 

In our Parallel Tempering simulations, several statisti- 
cally independent copies of the system at different en- 
ergies are simulated (fixed energies rather than fixed 
tem per atures, as it is normally performed in standard 

ptHhH). 

Each Monte Carlo time unit consists of two steps: 

1. For each copy of the system, we perform 10 5 x 
N particle move attempts at fixed energy (either 
single-particle displacements or particle-swap at- 
tempts). During this stage, each copy of the system 
is completely independent from the others. 

2. Copies of the system at neighboring energies try to 
exchange their particle configuration. We first try 
to sweep the two configurations at the lowest en- 
ergy, afterwards the second lowest with third low- 
est, etc. In this way, the particle-configuration at 
the lowest energy has a chance of getting to the 
highest energy in a single sweep. 

For the sake of clarity let us name A, B the two 
systems that are currently attempting to exchange 
their particle configuration. The exchange is ac- 



cepted with probability 
. [ ^({r^};eW)a, Af ({r| fl >};e^) l 

mm 1 , T-TT T^T . 1181 

L u N ({ri A) y,eW)u N ({rl B) y,e(B))\ 

The microcanonical weights ujn are given in 
Eq.®. 

Further details on the simulation are summarized in Ta- 
ble u 

Let us finally note that the here used Monte Carlo 
method is quite similar to that of Refs. [H, Hi!- We 
briefly mention the main differences. First, particle swap 
at fixed energy was not used in Refs. [12, [23] ■ Second, 
phase coexistence (and the related Maxwell construction) 
was not studied. Third, in the formulation of [27j . one 
has a single copy of the system that performs a random- 
walk in energy space: it is a sort of simulated annealing 
simulation [25| . rather than our parallel tempering. Be- 
sides, the approximation f3(e) rj (N — 2)/[2N((e — u))] 
is used, which coincides with Eq. only up to correc- 
tions of order 1/N. The formulation of [12| is somehow 
intermediate between simulated annealing and parallel 
tempering. The energy range of interest is spliced into 
non-overlapping subranges. Each copy of the system is 
assigned to an energy subrange, where it performs a sim- 
ulated annealing. From time to time one uses parallel 
tempering to exchange the copies of the system attached 
to neighboring energy subranges. 



C. Thermalization checks 

A crucial issue of PT simulations is to ensure thermal- 
ization. Fortunately, a nice feature of PT is that it pro- 
vides a sound thermalization check by controlling that all 
systems visit uniformly the whole range of energies [28| . 

At variance with the Potts case, phase coexistence in- 
side the energy gap between a fluid and a solid is apparent 
from the pdf of the quantity J 7 , defined in Eq. (fT4"f . Our 
results are shown in Fig. [TJtop). At values of e close 
to the transition, we identify two coexisting peaks. One 
of them is located at T ~ 1/-/V, as expected for an ho- 
mogenous fluid phase. On the other hand, the position 
of the large J- maximum becomes iV-independent (this is 
clearer at lower energies, see bottom panel in Fig. [I]), as 
it should occur for an inhomogeneous solid. Such phase 
coexistence makes us to expect a large growth with N 
of the autocorrelation times |29|. Actually, the pdf for T 
at low energies (Fig. [T] bottom) displays a shoulder at 
large J 7 , which corresponds to even more inhomogeneous 
solids. Hence, the PT dynamics is ruled by two differ- 
ent processes: tunneling from fluid to solid, and a second 
tunneling to even more inhomogeneous configurations. 

The random- walk in the energy space is best described 
through a PT time autocorrelation function (see Ref. (28[ 
for details) , that indeed can be fit to a double exponential 
for N = 256 and N = 500, see Fig.d Mind that the time 
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FIG. 1: pdf of J-, Ea. ()14[) at various representative values of 
e. Data in the top panel are computed at energy densities in 
the energy gap between the fluid and the solid phases. The 
double peak structure reveals phase coexistence (the position 
of the leftmost peak scales as 1/N). Data in the bottom panel 
are computed for e in the solid phase (the e-dependency there 
is very mild). 




FIG. 2: The (connected) time autocorrelation function for the 
energy in the PT can be fitted (dotted line) as (e(i') e(t'+t)} — 



a F se~ t/TFS +a S se~ t/TS ' 




1.05 



1.10 



1.15 



FIG. 3: Finite size effects in the Maxwell construction. Main 
panel: the inverse temperature /3(e) as a function of the en- 
ergy density e for various sizes of the sample. Inset: 5-/3 
phase diagram of the system obtained from the data in Ref.@. 



in this correlation functions correspond to the time-unit 
defined in Sect. IIII Bl It is not related to any physical 
time-correlation. 

As expected from the above discussion, we identify two 
different time scales in Table HI one associated to the co- 
existence of the homogeneous and inhomogeneous phase, 
Tpsy and a larger time, T55, related to the more inho- 
mogeneous configurations. For N — 864, we could only 
identify the tfs scale. Probably, tss is larger than the 
total time in our simulation. We remark that tfs for 
N — 256 can be estimated with a 5% accuracy, while only 
the order of magnitude of T55 is determined. We have 
explicitly checked that the effects of these very inhomo- 
geneous configurations on the Maxwell construction is 
fortunately smaller than our statistical errors 35| . Fur- 
thermore, from the point of view of our measured crys- 
talline order parameters (see below), the more inhomo- 
geneous configurations are not distinguishable from the 



main peak in the pdf. 



IV. NUMERICAL RESULTS 



The Maxwell construction 



As was mentioned in Sec. Ill B] in a microcanonical sim- 
ulation, a quantity of major interest is the (inverse) tem- 
perature, /3(e), see Eq. (jlll) . Thermodynamic stability 
requires that /3(e) be a decreasing function (i.e. positiv- 
ity of the specific heat). Yet, see main panel in Fig. [31 
this is not the case close to a first-order phase transition. 
The lack of monotonicity can be used to obtain the crit- 
ical temperature, surface tension, etc. through Maxwell 
construction (see below, and Ref.[20j for details). Gen- 
erally speaking, /3(e) has two distinct branches, one de- 
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scribing the fluid and the other the solid phase, where 
the specific heat C v = —(3 2 de/df3 is positive, connected 
by a thermodynamically instable line where C v < 0. Al- 
though at finite N the system does not undergo a real 
phase transition, there are various criteria to define an 
(inverse) critical temperature, ft , where the two differ- 
ent phases coexist with the same thermodynamic weight. 
Here we utilize the Maxwell construction, which amounts 
to obtain ft as a solution of: 



= 



r (/3 c N ) 



de (/3(e) - 



(19) 



where the energy e^{ft) (e%(ft)) in turn corresponds 
to the rightmost (leftmost) root of the equation /3(e) = 
ft. Eq. lfTBI shows that the Maxwell constructions 
amounts to the famous equal-height rule for the canonical 
probability-distribution function -Pg(e). 

In Fig. [3] we show the function /3(e) for N = 
256,500,864. At odds with other models displaying a 
first order transition, as N grows, both the supercooled 
fluid (fluid branch with /3 > ft) and the overheated solid 
(solid branch with /3 < ft) lines become longer. 

As for the values of ft reported in Table HI they 
decrease with N. Asymptotically, finite N corrections 
are of order 1/N (see [20] and references therein). A 
fit ft = /3 C °° + ai/N fails badly the \ 2 test. In other 
words, our estimates for ft are accurate enough to re- 
solve subleading scaling corrections in 1/N. Thus, we 
have used a different approach. Let us assume that scal- 
ing corrections take the form of a smooth function in 
1/N, ft = /3 C °° + a t /N + a 2 /N 2 + .... If we have at our 
disposal three values of A, we may compute a quadratic 
estimator (exact, up to corrections of order 1/N 3 ): 



A 



oo,quad 



A 



A'l 



+ A 



N 2 



+ A 



JV 3 





(Ax 


-N 2 )(N! 


-N 3 ) 




NI 




(A 2 


-Ni)(N 2 


— A3) 




Ni 




(A 3 


-Nx)(N 3 


-N 2 ) 



(20) 



Computing the statistical error in p^°-<i uad is trivial, since 
A^ 1 , ft 3 an d A^ 3 are statistically independent random 
variables. Using the data in Table Q] we get 



A 



00, quad 



4.624(20) , r ( 



00, quad 



1.4664(15). (21) 



However, the quadratic polynomial in 1/N that interpo- 
lates our values ft 1 , ft 2 and ft 3 displays a maximum 
by N w 256, and decreases for smaller N. Hence, / g^°^ uad 
probably overemphasizes curvature effects. On the other 
hand, a linear (in 1/N) extrapolation from Ni — 864 and 
N 2 = 500 yields 



A 



00, linear 



4.791(11), r ( 



00, linear 



1.4795(9). (22) 



The correct thermodynamic limit probably lies in be- 
tween of the two estimators r c °°' quad and r^°' lincal , above 
the kinetic glass transition at r g — 1.455(5). 



Furthermore, /3(e) also allows us to compute the sur- 
face tension [20(, 



= f N C de (/3(e) - ft ) , (23) 



[recall that equation /3(e) = ft has three solutions 

e^(AT ' ■ 

ble|H 



e%{ft) < e* N (ft) < e%(ft )}. Data is shown in Ta- 



B. Fractionation and crystalline ordering 

The need for generalized order parameters, Eq. (|16|) . 
follows from visual inspection of a typical N = 864 low- 
energy configuration, see Fig. 2) In fact, the smallest 400 
particles (particle index i < 400) and some of the inter- 
mediates (i € [600, 725]) show no sign of spatial order 
(bottom), while particles with i > 725 and i g [400,600] 
form crystalline planes. Ordered and disordered particles 
fill different regions of the sample. 

Our results in Fig. [5] confirm this picture. For x < 0.45 
the crystalline order parameters decay as 1/y/N, while 
for x — 0.55 and x — 0.95 we obtain results roughly N 
independent. Thus, while the latter group of particles 
form a crystal (Qe is somewhat smaller than expected 
for FCC ordering), the former one remains amorphous. 
As for polydispersities, in the two-components crystal we 
estimate that S ~ 0.15, while in the fluid 5 ~ 0.24. 



V. CONCLUSIONS 

In summary, we have studied in the microcanonical 
ensemble a soft-spheres model for liquids and colloids 
with a 24% polydispersity. Extrapolating by finite size 
scaling (FSS) to the thermodynamic limit the results ob- 
tained from the Maxwell construction in finite systems, 
we show that the critical temperature for the amorphous- 
crystal phase-separation is below the dynamic glass tran- 
sition, which makes dynamically difficult (although not 
impossible 10]) to observe such phase-separation. 

At low temperatures the system divides spatially into 
an amorphous and a crystalline part, in agreement with 
previous findings Q. The phase-separated amorphous is 
a stable fluid below its dynamic glass temperature, which 
is an optimal candidate to suffer a thermodynamic glass 
transition. On the other hand, the phase-separated solid 
displays crystalline order. Polydispersities on the coexist- 
ing amorphous and solid are smaller than in the fluid. In 
fact, particles distribute spatially according to their size 
following a complex pattern not described by any frac- 
tionation scenario known to us. However, we must men- 
tion that there are strong similarities with the results of 
very recent isobaric semigrand canonical simulations [30| . 
Although restricted to smaller system sizes (N = 256) 
and polydispersities (5 < 7% in the solid phase), these 
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n /3 C r c 


£ Pc <70 


TFS 


rss 




256 5.665(3) 1.5428(2) 
500 5.432(5) 1.5267(3) 
864 5.162(4) 1.5073(3) 


0.0035(2) 
0.0088(4) 


317(15) 
~1000 
~7000 


-20000 
-15000 


5x 32000 
2x30000 
1 x 12000 



TABLE I: Parameters of simulations and Maxwell construction. For each number of particles, N, we estimate two characteristic 
time scales tfs and tss for the PT random walk in energy (see text), in units of PT attempts. We perform 10 5 N MC steps at 
fixed energy, then try a PT sweep. We also report the total length of our simulations in units of PT sweeps (5 x 32000 stands 
for 5 independent runs of 32000 PT sweeps each). The energies chosen for the PT were evenly spaced e^+i — = 0.01, in the 
intervals [0.95, 1.14] (N = 256), [1.05, 1.2] (N = 500) and [1.08, 1.19] (N = 864). For N = 864 we added to the PT energy list 
the values 1.115, 1.125, 1.135 and 1.145 in the fluid-solid energy gap. We also report the (inverse) critical temperature (and the 
associated r c = as well as the dimensionless surface tension £/3 c v Oq, as computed from Maxwell's construction. 
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FIG. 4: Snapshot of a typical low energy configuration (TV = 864, e= 1.01). Top-left: whole system. Top-right: particles with 
index i > 725 and i £ [400, 600] . Bottom-left: particles i < 400. Bottom-right: particles i € [600, 725] . The size of the circles are 
proportional to the particle sizes. 



authors find as well that in the crystal phase the correla- 
tions between the fluctuating local particle-sizes extend 
to quite long spatial distances. 
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